# Replication package for 
# "The Economic Leverage of International Organizations in Interstate Disputes"
# Johannes Karreth
# June 30, 2017
# jkarreth@ursinus.edu

# This file: 14_claims_dy_analyses.R
# Purpose: Analyses of using force in claims, claim-year level

## Read in IGO & claim data

rm(list = ls())

# setwd("...")

library("rio")
library("dplyr")
library("rstanarm")
library("texreg")
library("xtable")
library("ggplot2")

# Source in functions
source("Functions/theme_jk.R")

dat <- import("claims_dyadyear.csv")
dat$host <- ifelse(dat$maxhost > 0, 1, 0)
dat$useforce <- ifelse(dat$maxhost > 3, 1, 0)
dat$war <- ifelse(dat$maxhost > 4, 1, 0)
dat$maxhost <- as.factor(dat$maxhost)
dat$cincdif <- log(abs(dat$cinc1 - dat$cinc2))
dat$atrade <- c(dat$expab / (dat$tottra + 0.00000001))
dat$btrade <-  c(dat$expba / (dat$tottrb  + 0.00000001))
dat$trademin <- apply(dat[, c("atrade", "btrade")], 1, min)
dat$trademax <- apply(dat[, c("atrade", "btrade")], 1, max)
dat$tradeimbal <- log(abs(dat$trademin - dat$trademax) + 0.0001)
dat$tradedep1 <- with(dat, {log((expab + expba) / cgdpca + 0.0001)})
dat$tradedep2 <- with(dat, {log((expab + expba) / cgdpcb + 0.0001)})
dat$tradedepmin <- apply(dat[, c("tradedep1", "tradedep2")], 1, min)
dat$tradedepmax <- apply(dat[, c("tradedep1", "tradedep2")], 1, max)
dat$gdppc_low_ln <- log(dat$gdppc_low)
dat$igo_ingr23_use = dat$igo_ingr2_use + dat$igo_ingr3_use

dat <- arrange(dat, claimdy, year)
dat <- mutate(group_by(dat, dyadid), 
              jointpol7_l1 = lag(jointpol7),
              atopally_l1 = lag(atopally),
              rivalry_th_l1 = lag(rivalry_th))

# All analyses only post-1945

dat <- dat[dat$year > 1945, ]

# Merge specific IGO data
igo_spec <- import("hl-igo_dyad.dta")

dat_spec <- merge(x = dat, y = igo_spec, by.x = c("dyadid", "year"), by.y = c("dyadid", "year"), all.x = TRUE, all.y = FALSE)
dat_spec$hligo_noWB <- dat_spec$igo_lev3_count_use - dat_spec$io2400_bin
dat_spec$hligo_noWB <- ifelse(dat_spec$hligo_noWB < 0, 0, dat_spec$hligo_noWB)
dat_spec$hligo_noIMF <- dat_spec$igo_lev3_count_use - dat_spec$io2880_bin
dat_spec$hligo_noIMF <- ifelse(dat_spec$hligo_noIMF < 0, 0, dat_spec$hligo_noIMF)
dat_spec$hligo_noWBIMF <- dat_spec$igo_lev3_count_use - dat_spec$io2880_bin - dat_spec$io2400_bin
dat_spec$hligo_noWBIMF <- ifelse(dat_spec$hligo_noWBIMF < 0, 0, dat_spec$hligo_noWBIMF)
dat_spec$hligo_noIFAD <- dat_spec$igo_lev3_count_use - dat_spec$io2760_bin
dat_spec$hligo_noIFAD <- ifelse(dat_spec$hligo_noIFAD < 0, 0, dat_spec$hligo_noIFAD)

# Analyses

# Outcome: Use of force

m1 <- Zelig::relogit(formula = useforce ~ igo_lev3_count_use_l1 + salint + saltan + terriss + jointpol7_l1 + rivalry_th_l1 + absidealdiff_l1 + atopally_l1, data = dat, bias.correct = TRUE, tau = mean(dat$useforce, na.rm = TRUE))
summary(m1)

m1_stanglmer <- stan_glmer(formula = useforce ~ igo_lev3_count_use_l1 + salint + saltan + terriss + jointpol7_l1 + rivalry_th_l1 + absidealdiff_l1 + atopally_l1 + (1 | claimdy), data = dat, family = binomial(link = "logit"), 
                           prior = normal(location = 0, scale = 10), 
                           prior_intercept = normal(0, 10),
                           seed = 20170207,
                           cores = 4)

summary(m1_stanglmer, pars = "beta")
saveRDS(m1_stanglmer, file = "Output_MCMC/claims_dy_m1.RDS")

# Outcome: At least one peaceful (bilateral and/or third party) attempt this year

m2 <- Zelig::relogit(formula = attanyp ~ igo_lev3_count_use_l1 + salint + saltan + terriss + jointpol7_l1 + rivalry_th_l1 + absidealdiff_l1 + atopally_l1, data = dat, bias.correct = TRUE, tau = mean(dat$attanyp, na.rm = TRUE))
summary(m2)

m2_stanglmer <- stan_glmer(formula = attanyp ~ igo_lev3_count_use_l1 + salint + saltan + terriss + jointpol7_l1 + rivalry_th_l1 + absidealdiff_l1 + atopally_l1 + (1 | claimdy), data = dat, family = binomial(link = "logit"), 
                           prior = normal(location = 0, scale = 10), 
                           prior_intercept = normal(0, 10),
                           seed = 20170207,
                           cores = 4)

summary(m2_stanglmer, pars = "beta")
saveRDS(m2_stanglmer, file = "Output_MCMC/claims_dy_m2.RDS")

# Tables

class(m1) <- "glm"
class(m2) <- "glm"
texreg(list(m1, m2),
        stars = 0.1,
        file = "Output_Tables-and-Figures/claims_dy_results.tex",
        reorder.coef = c(2:9, 1),
        custom.note = "* p < 0.05, one-tailed tests.",
        custom.model.names = c("Using force", "Peaceful attempts"),
        custom.coef.names = c("Intercept",
                              "IGOs with high leverage (lagged)",
                              "Intangible salience of claim", 
                              "Tangible salience of claim", 
                              "Territorial claim", 
                              "Joint democracy (lagged)", 
                              "Strategic rivalry (lagged)",
                              "UNGA ideal point difference (lagged)",
                              "Allies (lagged)"),
        digits = 3,
        caption = "Determinants of using force and peaceful bilateral settlement attempts in claim-years: logistic regression estimates (fit with maximum likelihood, using rare events bias correction)",
        caption.above = TRUE)

m1_tab <- m1_stanglmer$stan_summary[, c(1, 3)]
m1_tab <- data.frame(m1_tab[c(2:9, 1, 215), ])
m1_tab$varname <- c("IGOs with high leverage (lagged)",
                     "Intangible salience of claim", 
                     "Tangible salience of claim", 
                     "Territorial claim", 
                     "Joint democracy (lagged)", 
                     "Strategic rivalry (lagged)",
                     "UNGA ideal point difference (lagged)",
                     "Allies (lagged)",
                     "Intercept",
                     "Log-posterior density")
names(m1_tab) <- c("Mean", "Standard deviation", "Variable")
m1_tab <- m1_tab[, c(3, 1, 2)]

m1_xtab <- xtable(m1_tab, digits = 2)
print(m1_xtab, file = "Output_Tables-and-Figures/claims_dy_m1_table.tex", include.rownames = FALSE)

m2_tab <- m2_stanglmer$stan_summary[, c(1, 3)]
m2_tab <- data.frame(m2_tab[c(2:9, 1, 215), ])
m2_tab$varname <- c("IGOs with high leverage (lagged)",
                    "Intangible salience of claim", 
                    "Tangible salience of claim", 
                    "Territorial claim", 
                    "Joint democracy (lagged)", 
                    "Strategic rivalry (lagged)",
                    "UNGA ideal point difference (lagged)",
                    "Allies (lagged)",
                    "Intercept",
                    "Log-posterior density")
names(m2_tab) <- c("Mean", "Standard deviation", "Variable")
m2_tab <- m2_tab[, c(3, 1, 2)]

m2_xtab <- xtable(m2_tab, digits = 3)
print(m2_xtab, file = "Output_Tables-and-Figures/claims_dy_m2_table.tex", include.rownames = FALSE)

# First differences (using Zelig)
m1.zout <- Zelig::zelig(useforce ~ igo_lev3_count_use_l1 + salint + saltan + terriss + jointpol7_l1 + rivalry_th_l1 + absidealdiff_l1 + atopally_l1, 
               model = "relogit", 
               tau = mean(dat$useforce, na.rm = TRUE),
               case.control = c("prior"),
               bias.correct = TRUE,
               data = as.data.frame(dat))

x.high <- Zelig::setx(m1.zout, igo_lev3_count_use_l1 = quantile(model.frame(m1)$igo_lev3_count_use_l1, prob = 0.9))
x.low <- Zelig::setx(m1.zout, igo_lev3_count_use_l1 = quantile(model.frame(m1)$igo_lev3_count_use_l1, prob = 0.1))
s.out <- Zelig::sim(m1.zout, x = x.low, x1 = x.high)
summary(s.out)

m2.zout <- Zelig::zelig(attanyp ~ igo_lev3_count_use_l1 + salint + saltan + terriss + jointpol7_l1 + rivalry_th_l1 + absidealdiff_l1 + atopally_l1, 
                        model = "relogit", 
                        tau = mean(dat$attanyp, na.rm = TRUE),
                        case.control = c("prior"),
                        bias.correct = TRUE,
                        data = as.data.frame(dat))

x.high <- Zelig::setx(m2.zout, igo_lev3_count_use_l1 = quantile(model.frame(m2)$igo_lev3_count_use_l1, prob = 0.9))
x.low <- Zelig::setx(m2.zout, igo_lev3_count_use_l1 = quantile(model.frame(m2)$igo_lev3_count_use_l1, prob = 0.1))
s.out <- Zelig::sim(m2.zout, x = x.low, x1 = x.high)
summary(s.out)


m2.fd <- data.frame(fd = s.out$sim.out$x1$fd[[1]],
                    pos = ifelse(s.out$sim.out$x1$fd[[1]] > 0, 1, 0))
sum(m2.fd$pos) / length(m2.fd$pos)

p <- ggplot(data = m2.fd, aes(x = fd * 100)) + geom_density(fill = "gray", color = "gray") + geom_vline(xintercept = 0, color = "white") + theme_jk() + annotate("text", x = 3, y = 0.05, label = "95.8% > 0", color = "white") + xlab("Change in Pr(Peaceful settlement attempt)\n  comparing dyads with 6 vs. 1 joint memberships in IGOs with high leverage") + ylab("") + scale_y_continuous(breaks = NULL)

ggsave(p, file = "Output_Tables-and-Figures/claims_dy_m2_fd.pdf", width = 5.5, height = 3)